Basketball: how likely is it to score?

680 μs

11.3 ms

When playing basketball we can ask ourself: how likely is it to score given a position in the court?

We use data from NBA games from the season 2006 - 2007. We will consider all types of shots.

7.1 μs

xxxxxxxxxx
 
begin
using CSV
using DataFrames
using StatsPlots
​
season_2007 = CSV.read("./data/seasons/shots_season_2006_2007.csv");
​
#shots = vcat(season_2007, season_2008, season_2009);
shots = season_2007
​
shots[:distance] = sqrt.( shots.x .^ 2 +  shots.y .^ 2);
shots[:angle] = atan.( shots.y ./ shots.x );
filter!(x-> x.distance>1, shots)
end;

42 s

But how we interpret the data?

In the sketch below we show the drawing of a basketball court, its dimensions and how to interpret the data in the table.

19.6 μs

10 s

So, the x and y have their origin at the hoop, and we compute the distance from this point to where the shot was made. Also, we will compute the angle with respect to the x axis, showed as θ int he previous sketch.

20.5 μs

We now plot where the shots where made:

5.9 μs

22.3 s

We see that the shots are very uniformly distributed near the hoop, except for distances very near to the hoop, to see this better, we plot the histograms for each axis, x and y.

9.5 μs

But we are interested in the shots that were scored, so we filter now the shots made and plot the histogram of each axis.

4.9 μs

xxxxxxxxxx
 
shots_made = filter(x->x.result==1, shots);

378 ms

1.6 s

241 ms

If we plot a 3d plot of the count we obtain the plot wireplot shown below.

8.4 μs

xxxxxxxxxx
 
begin
    using StatsBase
    h = fit(Histogram, (shots_made.x, shots_made.y), nbins=40)
    plot(h) # same as histogram2d
    wireframe(midpoints(h.edges[2]), midpoints(h.edges[1]), h.weights, zlabel="counts", xlabel="y", ylabel="x", camera=(40,40))
title!("Histogram of shots scored")
end

2.8 s

The first model we are going to propose is a Bernoulli model.

Why a Bernoulli Distribution?

A Bernoulli Distribution results from a experiment where we have 2 possible outcomes, one that we usually called a success and another called a fail. In our case our sucess is scoring the shot and the other possible event is failing it.

The only parameter in the bernoulli distribution is the probability p of having a success. We are going to model this parameter as a logistic function:

21.2 μs

xxxxxxxxxx
 
begin
using Turing
using StatsFuns: logistic
end

13.5 s

1.9 s

Why a logistic function?

We are going to model the probability of shoot as a function of some variables, for example the distance to the hoop, and we want that our probability of scoring increases as we are getting closer to it. Also out probability needs to be between 0 an 1, so a nice function to map our values is the logistic function.

8.5 μs

So, the model we are going to propose is:

$p \sim l o g i s t i c (a + b * d i s t a n c e [i] + c * a n g l e [i])$

$o u t c o m e [i] \sim B e r n o u l l i (p)$

6.3 μs

But what values and prior distributions are we going to propose to a, b and c?

Let's see:

8.7 μs

Prior Predictive Checks: Part I

5 μs

Suppose we say that our prior distributions for a, b and c are going to be 3 gaussian distributions with mean 0 and variance 1. Lets sample and see what are the possible posterior distributions:

$a \sim N (0, 1)$

$b \sim N (0, 1)$

$c \sim N (0, 1)$

8.9 μs

201 ms

490 ms

We see that some of the predicted behaviours for p don't make sense. For example, for positive values of b, we are saying that as we increase our distance from the hoop, the probability of scoring also increase. So we propose instead the parameter b to be the negative values of a LogNormal distribution. The predicted values for p are shown below.

9.3 μs

So,m our model now have as prior:

$a \sim N o r m a l (0, 1)$

$b \sim L o g N o r m a l (1, 0.25)$

$c \sim N o r m a l (0, 1)$

7.5 μs

and sampling values for the above prior distribution and we obtain the plot shown below for the predicted values of p.

8.9 μs

xxxxxxxxxx
 
begin
b_prior_sampling_negative = rand(LogNormal(1,0.25), n_samples)
predicted_p_inproved = []
for i in 1:n_samples
        push!(predicted_p_inproved, logistic.(a_prior_sampling[i] .- b_prior_sampling_negative[i].*possible_distances))
end
    end

39.2 ms

105 ms

Now that we have the expected behaviour for p, we define our model and calculate the posterior distributions with our data points.

7.1 μs

Defining our model and computing posteriors

6.7 μs

Now we define our model to sample from it:

7.3 μs

logistic_regression (generic function with 1 method)

xxxxxxxxxx
 
@model logistic_regression(distances, angles, result,n) = begin
  N = length(distances)
  # Set priors.
    a ~ Normal(0,1) 
    b ~ LogNormal(1,0.25)
    c ~ Normal(0,1)
    for i in 1:n
        p = logistic( a -  b*distances[i] + c*angles[i])
        result[i] ~ Bernoulli(p)
    end
end

158 μs

1.4 μs

chain

Chains MCMC chain (1500×12×3 Array{Float64,3}):

Iterations        = 1:1500
Thinning interval = 1
Chains            = 1, 2, 3
Samples per chain = 1500
parameters        = a, b, c
internals         = acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, is_accept, log_density, lp, n_steps, nom_step_size, step_size

Summary Statistics
  parameters      mean       std   naive_se      mcse         ess      rhat  
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64  
                                                                             
           a    0.2233    0.1350     0.0020    0.0023   3648.1578    1.0000  
           b    1.7963    0.2827     0.0042    0.0045   4701.3859    1.0000  
           c   -0.0358    0.1205     0.0018    0.0062    362.7575    1.0080  

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%  
      Symbol   Float64   Float64   Float64   Float64   Float64  
                                                                
           a   -0.0113    0.1370    0.2236    0.3083    0.4751  
           b    1.2820    1.6064    1.7866    1.9737    2.3572  
           c   -0.1805   -0.0820   -0.0312    0.0160    0.1089

xxxxxxxxxx
 
# Sample using HMC.
chain = mapreduce(c -> sample(logistic_regression(shots.distance[1:n] ./ maximum(shots.distance[1:n] ), shots.angle[1:n], shots.result[1:n], n), HMC(0.05, 10), 1500),
    chainscat,
    1:3
)
​

80.3 s

xxxxxxxxxx
 
plot(chain, dpi=60)

4.6 s

xxxxxxxxxx
 
begin
    a_mean = mean(chain[:a])
    b_mean = mean(chain[:b])
    c_mean = mean(chain[:c])
end;

107 ms

Now plotting the posterior distribution for an angle of 45°, we have:

6 μs

xxxxxxxxxx
 
begin
p_constant_angle = []
    for i in 1:length(chain[:a])
        push!(p_constant_angle, logistic.(chain[:a][i] .- chain[:b][i].*possible_distances .+ chain[:c][i].*π/4));
        end
    end

591 ms

600 ms

The plot shows that the probability of scoring is higher as our distance to the hoop decrease, which makes sense.

5.3 μs

xxxxxxxxxx
 
begin
p_constant_distance = []
    for i in 1:length(chain[:a])
        push!(p_constant_distance, logistic.(chain[:a][i] .- chain[:b][i].*0.5 .+ chain[:c][i].*possible_angles));
        end
    end

596 ms

948 ms

We see that the model predict a almost constant probability for the angle.

4.5 μs

New model and prior predictive checks: Part II

3 μs

Now we propose another model with the form:

$p \sim l o g i s t i c (a + b^{d i s t a n c e [i]} + c * a n g l e [i])$

But for what values of b the model makes sense?

6.9 μs

xxxxxxxxxx
 
begin
    f1(x) = 0.3^x
    f2(x) = 1.5^x
    f3(x) = -0.3^x
    f4(x) = -1.5^x
end;

67 μs

1 s

Analysing the possible balues for b, the one that makes sense is f1, since we want increasing influence of the distance in the values of p as the distance decreases, since the logistic function has higher values for higher values of x.

9 μs

835 ms

Defining the new model and computing posteriors

5 μs

logistic_regression_exp (generic function with 1 method)

xxxxxxxxxx
 
@model logistic_regression_exp(distances, angles, result,n) = begin
  N = length(distances)
  # Set priors.
    a ~ Normal(0,1) 
    b ~ Beta(2,5)
    c ~ Normal(0,1)
    for i in 1:n
        p = logistic( a +  b .^ distances[i] + c*angles[i])
        result[i] ~ Bernoulli(p)
            end
end

97.6 μs

chain_exp

Chains MCMC chain (1500×12×3 Array{Float64,3}):

Iterations        = 1:1500
Thinning interval = 1
Chains            = 1, 2, 3
Samples per chain = 1500
parameters        = a, b, c
internals         = acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, is_accept, log_density, lp, n_steps, nom_step_size, step_size

Summary Statistics
  parameters      mean       std   naive_se      mcse        ess      rhat  
      Symbol   Float64   Float64    Float64   Float64    Float64   Float64  
                                                                            
           a   -0.9250    0.1393     0.0021    0.0085   218.1599    1.0281  
           b    0.1441    0.1068     0.0016    0.0065   216.2263    1.0316  
           c   -0.0333    0.0902     0.0013    0.0033   792.6944    1.0039  

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%  
      Symbol   Float64   Float64   Float64   Float64   Float64  
                                                                
           a   -1.1927   -1.0200   -0.9240   -0.8303   -0.6463  
           b    0.0082    0.0618    0.1214    0.2019    0.4080  
           c   -0.1804   -0.0793   -0.0298    0.0187    0.1105

xxxxxxxxxx
 
# Sample using HMC.
chain_exp = mapreduce(c -> sample(logistic_regression_exp(shots.distance[1:n] ./ maximum(shots.distance[1:n] ), shots.angle[1:n], shots.result[1:n], n), HMC(0.05, 10), 1500),
    chainscat,
    1:3
)
​

73.6 s

xxxxxxxxxx
 
plot(chain_exp, dpi=60)

22.6 ms

xxxxxxxxxx
 
begin
p_exp_constant_angle = []
    for i in 1:length(chain_exp[:a])
        push!(p_exp_constant_angle, logistic.(chain_exp[:a][i] .+ chain_exp[:b][i].^possible_distances .+ chain_exp[:c][i].*π/4));
        end
    end

741 ms

570 ms

555 ms

im3

3 s

Does the Period affect the probability of scoring?

4.2 μs

Now we will try to answer this question. We propose then a model, with parameter for each one of the four possible periods.

3.7 μs

logistic_regression_period (generic function with 1 method)

xxxxxxxxxx
 
@model logistic_regression_period(distances, angles, result, period, n) = begin
  N = length(distances)
  # Set priors.
    a ~ filldist(Normal(0,1), 4)
    b ~ filldist(Beta(2,5), 4)
    c ~ filldist(Normal(0,1), 4)
    for i in 1:n
        
        p = logistic( a[period[i]] +  b[period[i]] .^ distances[i] +                                 c[period[i]]*angles[i])
        
        result[i] ~ Bernoulli(p)
            end
end

116 μs

How the time left to end the period affects the probability of scoring?

3 μs